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Abstract 



We address the question of whether solids may be distinguished from fluids by their 
response to shear stress. 
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I. Introduction 

Our focus is the theoretical modeling, within statistical mechanics, of the 
solid/fluid phase transition of matter in thermal equilibrium, for instance the 
ice/water transition, at high pressure and temperature, and in particular the use of 
rigidity to distinguish the phases. 

There are no analytic proofs of a solid/fluid transition in any statistical me- 
chanics model which uses particles in space interacting through simple short range 
forces (see however [1]), though there are many simulations showing the transition, 
both Monte Carlo and molecular dynamics. Since we concentrate on the transi- 
tion at high pressure and temperature, at which short range repulsion dominates 
the interparticle interaction, the classical hard sphere model is appropriate [2, page 
84]; again there are no proofs of a phase transition in this model, but there are 
many simulations [2,3]. Traditionally such a transition is "understood" theoreti- 
cally through an order parameter associated with some global (emergent) property 
of the bulk material, in particular the molecular-level crystalline symmetry which 
can be detected experimentally in X-ray scattering [4] . This paper follows an alter- 
native proposal of Anderson [5], namely the use of rigidity, the response to stress, 
to distinguish the phases, for instance in a hard sphere model. 

Stress, both pressure and shear, will be understood here as an external influ- 
ence (force per unit area) on the boundary of a finite sample of the material, with 
pressure acting on the volume and shear on the shape; its extension to a uniform 
stress field inside the material is an important issue to be addressed. Pressure is 
commonly modeled in statistical mechanics as a parameter in a pressure ensemble 
which controls the deformations of volume. A similar approach can be followed for 
shear stress, using a set of parameters controlling fluctuations of the shape of the 
container, though this is much less common; see however [6,7]. 

One might in principle be able to model the distinction between ice and water 
through statistical mechanics estimates of compressibility (— [dV/dp]/V), in which 
V is volume and p is pressure. Unfortunately the compressibility of ice and water 
are both high and the difference is relatively small, a common circumstance for a 
solid/fluid transition. However the difference in the corresponding elastic shear con- 
stants is, experimentally, much greater, since elastic shear constants are negligible 
for (isotropic) fluids. This suggests an advantage in using shear instead of pressure 
to distinguish a solid from its melt. 

But, as emphasized for instance in the recent paper of Sausset et al [8], for a 
material in equilibrium any linear response to macroscopic shear must be transient 
in time, making it harder to model an elastic shear constant within equilibrium sta- 
tistical mechanics. Indeed, there are proofs that in equilibrium statistical mechanics 
the free energy is independent of the shape of the material [9,10], which suggests 
that shear stress be properly considered as taking a material out of equilibrium. In 
[8] solids are treated as highly viscous fluids and solids are distinguished from fluids 
by a dynamical feature, the divergence of the viscosity as the shear stress rate van- 
ishes. Since we are using an equilibrium model we focus instead on a spatial issue, 
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namely the question of whether the response to shear is localized at a small part of 
the material or is (uniformly) distributed throughout the material. We concentrate 
on the response to shear stress in the limit of zero shear, computed before any bulk 
limit. (At zero shear stress the shape of the material is unconstrained; imagine a 
triaxial shear cell, with negligible friction and in zero gravity.) Mathematically, we 
are taking advantage of the fact that the limit of vanishing shear need not commute 
with the infinite volume limit. 

Our approach is based on the following idea. A configuration of hard spheres 
at high pressure must be approximately crystalline, with most particles trapped in 
a cage of neighbors. A macroscopic change in container shape can be accommo- 
dated by adjustments only near the boundary, without affecting the structure in 
the bulk interior; see Fig. 1. On the contrary a very small change of shape cannot 
be accommodated in this way and might result in a small deformation of internal 
structure throughout the configuration rather than just near the boundary; see Fig. 
2. 
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Figure 1. The result of a macroscopic change of strain: the interiors of A) and B) 
are the same. 
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Figure 2. The result of a small change of strain: the interiors of A) and B) are 
different, the underlying lattice becoming distorted. 

In other words, there may be a regime in which the response to shear is linear 
and extends throughout the material, but such a regime, measured by the angle 
of deformation, would have to vanish with the size of the system, constituting an 
equilibrium alternative to the dynamical effect discussed by Sausset et al. If indeed 
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the response extends throughout the material, which is by no means evident, it 
would be appropriate to measure it, in a finite system at constant high pressure, by 
the rate of change in density <p with respect to shear stress /, computed at / = 0. 
We expect this to be large in magnitude at high pressure, in the solid phase. At 
low pressure the model should represent a fluid with negligible response, and we 
propose this difference as a means of detecting the solid/fluid phase transition. 

In the remainder of this paper we give indirect support to the proposition 
that for infinitesimal shear the material responds linearly throughout the finite 
sample. We do this by simulation of the two dimensional hard disk model, in 
a stress (pressure and shear) ensemble. Our measurements in this model show 
an emerging resistance to shear at volume fraction about 0.7, very close to the 
known transition(s) for hard disks. Our simulation does not show the details of 
the transition, which are well displayed in the recent tour-de-force by Bernard and 
Krauth [11]. Instead, the point of this work is merely to illustrate the feasibility of 
using density response to infinitesimal shear to probe a solid/fluid phase transition, 
in the tradition of Anderson [5] . 

II. The Model and Simulations 

We consider arrangements of a fixed number N of hard disks of fixed radius 
a inside various parallelograms, with the volume and shape of the parallelograms 
allowed to vary. More precisely, we consider arrangements of such disks inside 
boundaries formed by placing disks along the edges of parallelograms at regular 
intervals as in Fig. 3. These boundary disks lie on a regular triangular lattice 
when the underlying parallelogram is rectangular, and all the parallelograms are 
related to each other via maps (on the boundary disk centers) of the form (x, y) — > 
(Xx + v\y, Xy) for real A, with A > 0. 

We employ a "stress ensemble" which uses parameters p and / to control the 
volume V and deformations of shape a, respectively, of the parallelograms. (We 
want to use the Lagrange multipliers p and / to model pressure and shear stress.) 
More precisely, we consider probability measures (states) which minimize the free 
energy 

F(pJ) = S - fipV + 0faV. (1) 

Here (3 is the inverse temperature and a is the angle of inclination of a parallelogram, 
with a = representing a rectangle (see Fig. 3). Such an ensemble has partition 
function oo oo / 

Z pJ = [ [ ( [ dC) exp(-(3pV + (3faV)dVda, (2) 

JO JO \Jv,a J 

where J v dC represents integration over all arrangements of hard disks in a par- 
allelogram of volume V and shape a. (Temperature plays a simple role in hard core 
models such as this, so we will assume that numerically (3 = 1 where convenient. 
Also we are using the usual "reduced" ensemble in which the velocity variables have 
been integrated out.) By the change of coordinates 

V>v> : (x, y) (V-^ix - tan(a)y), V~ 1 / 2 y), (3) 
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considered as a function on the disk centers, the partition function may be rewritten 

/ / / ^vl(Q)) vN eM-PpV + PfaV)dQdVda, (4) 
Jo Jo Jq 



as 

Z pJ 



where $(C) = 1 if no two disks of radius a with centers from C overlap, $(C) = 
otherwise, and Q C M 2 is a fixed rectangular area. The probability density of an 
arrangement of ( no nover lapping) hard disks in a parallelogram of volume V and 
shape a is then proportional to 



V" exp{~PpV + pfaV). (5) 

Because we are interested only in infinitesimal shear, we impose the restriction 
< a < 0.01, with a measured in radians. (For the densities and a we consider, 
the boundary disks do not come close to overlapping.) 




Figure 3. An arrangement of disks in a parallelogram. Boundary disks are in bold; 
a is the angle formed between the boundary disks and a vertical line. 

Let 4> p j be the average volume fraction of arrangements at fixed p and /. To 
measure the volume response of arrangements of disks to an infinitesimal change in 
shape, we estimate the derivative 



r(p) := 



(6) 

f=o 



Of 

By definition the average volume fraction (p p j is given by 

,j= / / ^y 1 a (Q))(Nna 2 /V)V N exp(-^pV + PfaV)dQdVda. (7) 
Jo Jo Jq 
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Differentiating with respect to /, one obtains 



r(p) = (5Nna 2 [{a) pfi - (i\ (7a)„, 



(8) 



with (-)pj representing an average value at fixed p and /. Applying equation (8), 
we obtain T(p) from the average values of a, 1/V and Va, in our simulations at 
pressure p with / = 0. 
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Figure 4. Time series for systems of iV = 10656 disks, with t in units of 3.5 x 10 8 
Monte Carlo steps, and at / = 0, clockwise from top left: a,b) Volume fraction <fi 
vs. t; c,d) shape a vs. t. 

The change from T(p) = to T(p) << represents the development of rigidity 
in the system, in the sense that the system has to overcome pressure p to expand 
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in response to infinitesimal shear. 

Our simulations use Monte Carlo steps consisting of moves which change the 
size and shape of the arrangements of disks, as well as so-called "event-chain" 
movements of multiple disks [11]. In the former types of moves, the coordinates 
of disk centers (x,y) are mapped to (Xx,Xy) or (x + py,y), respectively. Here 
A = (V + rje) 1 / 2 /V 1 / 2 and p = u8, where 77, v are (independent) random variables 
distributed uniformly in (—1, 1), and e and 5 are positive real parameters. 




Figure 5. Volume fraction data, clockwise from top left: a) Average volume frac- 
tion <p Py o vs. pressure (3p(2a) 2 , with confidence intervals smaller than the data 
markers; b) Differential volume response T vs. pressure (3p(2a) 2 ; c) Differential 
volume response F vs. average volume fraction Pj o; d) Mixing time, as a fraction 
of total number of Monte Carlo steps, vs. P; o- 
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If such a move results in overlap of disks, it is rejected. These types of moves are 
(each) attempted with frequency N~ 1 / 2 /A. In the latter type of move, employed 
recently in [11], a non-boundary disk and a random direction are selected, with the 
direction being up, down, left or right (that is, parallel to one of the coordinate 
axes). Additionally a displacement L is selected uniformly at random from the 
interval (0, yV — Nrra 2 /2). The particle is then moved along the chosen direction 




0.005 



-0.005 



- U -0.01 



8.5 9 9.5 

P P (2a) 2 




-0.015 



-0.02 



-0.025 



0.01 



0.005 



-0.005 



-0.01 



-0.015 



-0.02 



-0.025 



-0.03 




P P (2°) 



Figure 6. Confidence intervals for differential volume response V vs. pressure 
/3p(2cr) 2 , for systems with (clockwise from top left) N = 3520, N = 5451, N = 7790, 
and N = 10656. 
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until it strikes another particle, at which point that particle moves in the same 
direction until it strikes another particle, and so on. The process continues until 
a total displacement of L is obtained, the total displacement being the sum of the 
displacements of all the particles. If the process results in the displacement of a 
boundary particle, then the move is rejected. (It is in principle possible that such 
moves result in disks moving outside the boundary, but this does not occur for the 
pressures we simulate.) Such moves are attempted with frequency 1 — N~ 1 / 2 /2. 

We investigate systems with N = 3520, N = 5451, N = 7790, N = 10656, and 
N = 13970 disks, beginning with perfectly crystalline arrangements of the disks. At 
each pressure p we run these systems to 2 x 10 10 , 2.5 x 10 10 , 3 x 10 10 , 3.5 x 10 10 , and 
4 x 10 10 Monte Carlo steps, respectively. This results in about 2 x 10 7 displacements 
per particle, and about 10 7 fluctuations in volume and shape, for each p and system 
size N. We checked that our runs were long enough for volume fraction <p and 
shape a to equilibrate, after a burn-in time of at most about 10% of the run (with 
the exception of the shape a at large p, as discussed below); see Fig. 4. Therefore 
in our main data, shown in Fig. 5 (with confidence intervals in Fig. 6), we have 
thrown away the first 10% of each Monte Carlo run. Along with our main data we 
also measured mixing times, defined as the number of Monte Carlo steps before the 
standard (unbiased) autocorrelation of the time series for </> P; o first crosses zero; see 
Fig. 5. Excluding the largest system (N = 13970), mixing times were no more than 
15% of our Monte Carlo runs. For 90% confidence intervals, we run 10 independent 
copies of every simulation and use Student's t-distribution with 9 degrees of freedom 
on the average values obtained from each of the 10 copies; see Fig. 6. We do not 
compute confidence intervals for the largest system, N = 13970. 

We find the volume response parameter T(p) defined in (8) exhibits the fol- 
lowing behavior (see Fig. 5). At low pressure p or volume fraction = 4> Pj q, Y(p) 
is approximately zero, indicating there is no volume response to an infinitesimal 
shear. We interpret this as meaning the hard disk fluid does not resist a small shear 
stress. As <fi rises above 0.70, the volume response V begins to deviate from zero. 
Our data is not fine enough to distinguish the details of the transitions shown in 
[11], and in particular does not justify estimating specific transition values for p. 

We note that the volume response Y measured in our simulations tapers off to 
zero at large p (or high density). We do not expect this tapering to be accurate; 
instead we interpret this as a sign that the simulations begin to get "stuck" as 
densities increase. This is confirmed by checking that the simulations no longer 
explore the full range of shapes a; see Fig. 4c) -d). We expect the true behavior of 
r, as a function of <fi, to be non-increasing, indicating a volume response into the 
nonfluid phases. We conclude, then, that the hard disk solid resists a small positive 
shear stress, while the fluid does not. 

III. Summary 

The rigidity of solids can be modeled in various ways. We have chosen to use 
equilibrium statistical mechanics, in large extent because it is the most convincing 
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formalism in which to distinguish solids from fluids, which is our motivation for 
studying rigidity, following Anderson [5]. Even within equilibrium statistical me- 
chanics one could model response to shear more simply with a harmonic crystal 
model [12: Chapter 22], with long range quadratic forces between particles assigned 
neighboring equilibrium sites. In fact this is commonly used to model sound (pres- 
sure) propagation, but does not exhibit a fluid phase and therefore does not allow 
comparison between solid and fluid phases, which is the purpose of our work. 

The most awkward consequence of using equilibrium statistical mechanics is 
that to obtain the sharp solid/fluid distinction one must take the infinite volume 
limit while for infinite systems one cannot model shear stress, as noted in the intro- 
duction. Our solution to this was to model the shear in finite volume - where there 
are no well defined solid or fluid phases but the system can model shear stress - 
and measure the volume response in the limit of vanishing shear, before taking the 
infinite volume limit. 

The other weakness of our approach is technical: in order to measure the 
volume response to stress we employed variation in both strain and volume, which 
is costly in simulation time compared to the usual Monte Carlo technique for the 
hard disk system, which uses fixed volume, strain and particle number. We favor 
this ensemble for its theoretical advantages: calculating response to stress in our 
ensemble is far simpler than in an ensemble which fixes, say, density and strain - 
in our ensemble we can calculate the density response directly from fluctuations, 
whereas an equivalent analysis in the latter model would require taking a numerical 
derivative of average pressure as shear strain vanishes, with pressure computed 
approximately (perhaps from a virial expansion; see [11]). Our data, due to the 
large computation time associated with fluctuations in volume and strain, is not 
sufficient to demonstrate the details of the transitions, as is done in [11]. We feel 
this is acceptable in exchange for demonstrating the feasibility of shear response as 
a theoretical tool to analyze the solid/fluid transition. 

In conclusion we note that our approach is similar to the analysis of the di- 
latancy transition recently found experimentally [13] in (nonequilibrium) static, 
granular matter, and its modeling [14] with a stress ensemble. In effect we are 
proposing to model the solid/fluid transition of equilibrium matter as a dilatancy 
transition, a sharp change between the solid and fluid equilibrium phases in their 
volume response to (infinitesimal) shear, instead of by a change in symmetry of the 
molecular pattern, as is the common practice. 

Acknowledgements. It is a pleasure to acknowledge useful discussions with Giulio 
Biroli, Daan Frenkel and Dan Stein. 
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